#load pakages
if(!require(pacman)){install.packages("pacman")}
pacman::p_load(tidyverse, haven, ggpubr, rstatix, glue, estimatr, ggplot2, MetBrewer)

#load data
data <- read_dta("data/hgopy_pilot_public.dta")
#create the new price variable with discounts grouped and counselling style as factors
data <- data %>% mutate(data, 
                        pdx = ifelse(price_larc_cat==3,'Full','Discounted'),
                        pdn = ifelse(price_larc_cat==3,0,1), 
                        pd = factor(pdx),
                        pd = fct_reorder(pd,pdn),
                        csx = ifelse(view==0,'MiM',
                              ifelse(view==1,'IDM',
                              ifelse(view==2,'SDM','None'))),
                        cs = factor(csx),
                        cs = fct_reorder(cs,view))

###Plot prep

#tests for the within group comparisons, Only MiM
tests_wncs_mim <- data %>%
  select(cs,pd,adopted_larc) %>%
  filter(cs=='MiM') %>% 
  group_by(cs) %>% 
  t_test(adopted_larc ~ pd, detailed = TRUE)
#gotta add the y position...
tests_wncs_mim <- tests_wncs_mim %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values and estimates nicely 
tests_wncs_mim <- tests_wncs_mim %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs_mim)

#tests for the within group comparisons, excluding MiM
tests_wncs <- data %>%
  select(cs,pd,adopted_larc) %>%
  filter(cs!='MiM') %>% 
  group_by(cs) %>% 
  t_test(adopted_larc ~ pd, detailed = TRUE)
#gotta add the y position...
tests_wncs <- tests_wncs %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values nicely 
tests_wncs <- tests_wncs %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs)

#tests for the between group comparisons, but excluding MiM
tests_bwcs <- data %>%
  select(cs,pd,adopted_larc) %>%
  filter(cs!='MiM') %>% 
  group_by(pd) %>%
  t_test(adopted_larc ~ cs, detailed = TRUE)
#gotta add the y position...
tests_bwcs <- tests_bwcs %>% 
  add_xy_position(x = "cs" , group = "pd" , dodge = 0.8) #fun = "mean_sd" ,
#format the p values nicely 
tests_bwcs <- tests_bwcs %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_bwcs)

#tests for the diffs-in-diffs
tests_did_lm <- data %>% 
  select(cs,pd,adopted_larc) %>%
  filter(cs!='MiM') %>% 
  lm_robust(adopted_larc ~ pd * cs , . , se_type = "stata")
print(summary(tests_did_lm))
didb <- summary(tests_did_lm)$coefficients["pdDiscounted:csSDM", "Estimate"] %>% 
  round(., digits=3)
didp <- summary(tests_did_lm)$coefficients["pdDiscounted:csSDM", "Pr(>|t|)"] %>% 
  round(., digits=3)
#might have to make a separate object with x y and the values we want
did <- data.frame(a = paste("DiD = ",didb," ","p = ",didp))
print(did)





###Plot

ggbarplot(data, x = "cs" , y = "adopted_larc", 
          fill = "pd", 
          add = c("mean_ci"),
          position = position_dodge(0.8),
          xlab = "Counselling style" , 
          title = "Fraction of clients who adopted a LARC",
          ylab = "",
          label = TRUE,
          lab.nb.digits = 3,
          lab.hjust = -.35,
          lab.size = 3) +
  #add in the p-values and differences, all default set at 1 so have to manually lower them down
  stat_pvalue_manual(tests_wncs, label = "d = {est}; p = {prounded}",
                     tip.length = 0.01,
                     bracket.nudge.y = -0.45,
                     label.size = 3) +
  stat_pvalue_manual(tests_wncs_mim, label = "d = {est}; p = {prounded}",
                     tip.length = 0.01,
                     bracket.nudge.y = -0.225,
                     label.size = 3) +
  stat_pvalue_manual(tests_bwcs, 
                     label = "d = {est}; p = {prounded}",
                     color = "black", 
                     tip.length = 0.01,
                     step.increase = 0.07,
                     bracket.nudge.y = -0.375,
                     label.size = 3) +
  #add the DID estimates 
  geom_text(data = did, 
            aes(x = 2.5,  y = .795, label = a),
  					size = 3) +
  annotate("segment", x = 1.8, xend = 3.2, y = .775, yend = .775, size=0.25) +
  annotate("segment", x = 1.8, xend = 1.8, y = .775, yend = .763, size=0.25) +
  annotate("segment", x = 3.2, xend = 3.2, y = .775, yend = .763, size=0.25) +
  #make it look nice
  theme_classic() +
  scale_fill_manual(values = met.brewer("Signac", 8)[c(3,7)]) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12),
        legend.text = element_text(size=12),
        legend.position= "bottom",
        legend.title = element_text(size=12),
        panel.grid.minor = element_blank(), 
        panel.grid.major = element_blank()) +
  guides(fill=guide_legend(title="LARC prices:")) 

  
ggsave("output/figures/Fig2-view-prices.png", dpi = 320, scale = 0.8, height = 14, width = 20, units = "cm")


 
